% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replication "Deconstructing the Yield Curve"
% Crump and Gospodinov (2024)
% Date: 26-JUL-2024
% Function for estimating VAR(1) model
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
function var1Out = var1BC(Y,B,biasCorrectExplosive)

[T,k] = size(Y);
y = Y(2:end,:);
X = [ones(size(y,1),1), Y(1:end-1,:)];
PsiOLS = (X'*X)\(X'*y);
vOLS = y - X*PsiOLS;
PhiOLS = PsiOLS(2:end,:)';

if (B > 0) && (max(abs(eig(PhiOLS))<1) || biasCorrectExplosive)
    bFlag = 1;
    bsPsiOLS = NaN([size(PsiOLS),B]);    
    for b = 1:B           
        idxB = randi(T-1,[T-1,1]);  
        tmpYB = buildVAR(PsiOLS(1,:)', PsiOLS(2:end,:)', vOLS(idxB,:)', Y(1,:)')';
        YB = [Y(1,:); tmpYB];
        bsy = YB(2:end,:);
        bsX = [ones(size(bsy,1),1), YB(1:end-1,:)];        
        bsPsiOLS(:,:,b) = (bsX'*bsX)\(bsX'*bsy);        
    end
    
    bsPsiOLS_mean = mean(bsPsiOLS,3);
    bsPsiOLS_diff = bsPsiOLS_mean-PsiOLS;

    flag = 1;
    kappa = 1;
    biasAdj = bsPsiOLS_diff;
    while (flag == 1)
       
        PsiBC = PsiOLS - biasAdj;

        if any(abs(eig(PsiBC(2:end,:)))>=1)
            kappa = kappa - .01;
            biasAdj = kappa*biasAdj;
        else
            flag = 0;
        end
    end

else
    PsiBC = PsiOLS;
    bFlag = 0;
end

var1Out.PsiOLS = PsiOLS;

var1Out.muOLS = PsiOLS(1,:)';
var1Out.PhiOLS = PsiOLS(2:end,:)';

var1Out.VOLS = [NaN(1,k); y - X*PsiOLS];

var1Out.SigOLS = (var1Out.VOLS(2:end,:)'*var1Out.VOLS(2:end,:))/(T-numel(PsiBC));

var1Out.PhiBC = PsiBC(2:end,:)';
if bFlag
    var1Out.muBC = mean(y)' - var1Out.PhiBC*(mean(X(:,2:end))');
else
    var1Out.muBC = var1Out.muOLS;
end
var1Out.PsiBC = [var1Out.muBC'; var1Out.PhiBC'];

var1Out.VBC = [NaN(1,k); y - X*var1Out.PsiBC];
tmp = var1Out.VBC - ones(T,1)*mean(var1Out.VBC, 'omitnan');

var1Out.SigBC = (tmp(2:end,:)'*tmp(2:end,:))/(T-numel(PsiBC));

if max(abs(eig(var1Out.PhiBC))>=1), var1Out.flag = 1; else, var1Out.flag = 0; end

